Download the files with R:
baseurl <- "http://earth-info.nga.mil/GandG/wgs84/gravitymod/egm2008/GIS/world_geoid"
fs <- file.path(baseurl, c("n45w180.zip", "n45w135.zip", "n45w90.zip", "n45w45.zip", "n45e00.zip", "n45e45.zip",
"n45e90.zip", "n45e135.zip", "n00w180.zip", "n00w135.zip", "n00w90.zip",
"n00w45.zip", "n00e00.zip", "n00e45.zip", "n00e90.zip", "n00e135.zip",
"s45w180.zip", "s45w135.zip", "s45w90.zip", "s45w45.zip", "s45e00.zip",
"s45e45.zip", "s45e90.zip", "s45e135.zip", "s90w180.zip", "s90w135.zip",
"s90w90.zip", "s90w45.zip", "s90e00.zip", "s90e45.zip", "s90e90.zip",
"s90e135.zip"))
for (i in seq_along(fs)) download.file(fs[i], basename(fs)[i], mode = "wb")
## unzip
for (i in seq_along(fs)) utils::unzip(basename(fs[i]))
## build the next string with R too
txt <- gsub(".zip", "", basename(fs))
cat(paste(file.path(txt, txt, "w001001.adf"), collapse = " "), "\n")
Mosaic them all with GDAL, and do a quick decimation, choose your own percentage. (This is longlat data so we really need more sophisticated resampling to get data evenly distributed on the globe, doable in R but for another time.)
gdalbuildvrt geoid.vrt n45w180/n45w180/w001001.adf n45w135/n45w135/w001001.adf n45w90/n45w90/w001001.adf n45w45/n45w45/w001001.adf n45e00/n45e00/w001001.adf n45e45/n45e45/w001001.adf n45e90/n45e90/w001001.adf n45e135/n45e135/w001001.adf n00w180/n00w180/w001001.adf n00w135/n00w135/w001001.adf n00w90/n00w90/w001001.adf n00w45/n00w45/w001001.adf n00e00/n00e00/w001001.adf n00e45/n00e45/w001001.adf n00e90/n00e90/w001001.adf n00e135/n00e135/w001001.adf s45w180/s45w180/w001001.adf s45w135/s45w135/w001001.adf s45w90/s45w90/w001001.adf s45w45/s45w45/w001001.adf s45e00/s45e00/w001001.adf s45e45/s45e45/w001001.adf s45e90/s45e90/w001001.adf s45e135/s45e135/w001001.adf s90w180/s90w180/w001001.adf s90w135/s90w135/w001001.adf s90w90/s90w90/w001001.adf s90w45/s90w45/w001001.adf s90e00/s90e00/w001001.adf s90e45/s90e45/w001001.adf s90e90/s90e90/w001001.adf s90e135/s90e135/w001001.adf
gdal_translate geoid.vrt geoid04.tif -outsize 4% 4% -co COMPRESS=LZW -co TILED=YES
Now back to R, read that simplified raster and convert to 3D.
library(rglwidget)
library(raster);library(rgdal)
## Loading required package: sp
## rgdal: version: 1.1-4, (SVN revision 594)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
## Path to GDAL shared files: E:/inst/R/R/library/rgdal/gdal
## GDAL does not use iconv for recoding strings.
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: E:/inst/R/R/library/rgdal/proj
## Linking to sp version: 1.2-1
r <- raster(file.path(gpath, "geoid04.tif"))
extent(r) <- extent(-180, 180, -90, 90)
library(gris) ## devtools::install_github("mdsumner/gris")
# Build geometry
geoid0 <- bgl(r, z = r * 5000)
geoid <- geoid0
## convert long-lat-z to XYZ
z <- t(llh2xyz(t(geoid$vb[1:3, ])))
## build map
library(rworldmap)
## ### Welcome to rworldmap ###
## For a short introduction type : vignette('rworldmap')
data(coastsCoarse)
globe <- gris(coastsCoarse)
## get natural earth image
u <- "http://naturalearth.springercarto.com/ne3_data/8192/textures/2_no_clouds_8k.jpg"
if (!file.exists(basename(u))) download.file(u, basename(u), mode = "wb")
tmp <- setExtent(brick(basename(u)), extent(-180, 180, -90, 90))
projection(tmp) <- "+proj=longlat +ellps=WGS84"
tmp <- projectRaster(tmp, raster(extent(tmp), nrow = nrow(tmp)/8, ncol = ncol(tmp)/8, crs = projection(tmp)))
imname <- gsub("jpg$", "png", basename(u))
if (!file.exists(imname)) writeGDAL(as(tmp, "SpatialGridDataFrame"), imname, drivername = "PNG")
And plot
library(colorspace)
##
## Attaching package: 'colorspace'
## The following object is masked from 'package:raster':
##
## RGB
pal <- colorRampPalette(c("blue", "lightblue", "white", "orange", "red"))(24)
scl <- function(x) (x - min(x, na.rm = TRUE)) / diff(range(x, na.rm = TRUE))
library(rgl)
##
## Attaching package: 'rgl'
## The following object is masked from 'package:gris':
##
## triangulate
shade3d(geoid, col = pal[scl(geoid0$vb[3,geoid$ib]) * (length(pal) - 1) + 1])
# for (i in seq(nrow(globe$o))) {
# xx <- mkpslg(globe[i, ])
# xx$XYZ <- llh2xyz(cbind(xx$P, extract(r, xx$P) + 1000 ))
# lines3d(xx$XYZ[xx$S, ])
# }
subid <- currentSubscene3d()
rglwidget(elementId="geoid")